#include <stdio.h>
#include <math.h>
#include <time.h>
#include "stdafx.h"
#include "cv.h"
#include "highgui.h"
#include "cxcore.h"

void decide_scale(double scales[3], const int n, int size)
{
	for(int i=0; i<n; i++)
		scales[i] = (double) 2.+(i)*(size/n);
}

void generate_coef3(double c[5], double sigma)
{
	double q;
	

	//decide q
	if(sigma>=2.5)
		q = 0.98711*sigma - 0.96330;
	else if(sigma<2.5 && sigma>=0.5)
		q = 3.97156 - 4.14554*sqrt(1-0.26891*sigma);
	else
		q = 0.1147705018520355224609375;

	double q2 = q*q;
	double q3 = q*q*q;

	//decide coefficients of recursive filter
	//b0, b1, b2, b3
	c[0] = 1.57825+(2.44413*q)+(1.4281*q2)+(0.422205*q3);
	c[1] = (2.44413*q)+(2.85619*q2)+(1.26661*q3);
	c[2] = -((1.4281*q2)+(1.26661*q3));						
	c[3] = 0.422205*q3;
	c[4] = 1 - ((c[1]+c[2]+c[3])/c[0]); //B
}


void gaussian_conv(double b[5], CvMat* in, CvMat* out)
{
	double invb0 = 1./b[0];
	int bufsize = in->width+3;
	std::vector<double> w1(bufsize);
	std::vector<double> w2(3*bufsize);

	CvScalar elem, elem1;
	for(int y=0; y<in->height; y++)
	{
		for(int i=0; i<3; i++)
		{
			elem = cvGet2D(in,y,0);
			//forward
			w1[0] = elem.val[i];
			w1[1] = elem.val[i];
			w1[2] = elem.val[i];

			for(int x=0; x<in->width; x++)
			{
				elem = cvGet2D(in,y,x);
				w1[x+3] = b[4]*elem.val[i] + ((b[1]*w1[x+2]+b[2]*w1[x+1]+b[3]*w1[x])*invb0);
			}

			//backward
			w2[i*(in->width+2)+in->width+0]=w1[in->width+2];
			w2[i*(in->width+2)+in->width+1]=w1[in->width+2];
			w2[i*(in->width+2)+in->width+2]=w1[in->width+2];
			for(int x=in->width-1; x>=0; x--)
			{
				w2[i*(in->width+2)+x] = b[4]*w1[x] + ((b[1]*w2[x+1]+b[2]*w2[x+2]+b[3]*w2[x+3])*invb0);
				//cvSetReal2D(out, x, y, w2[x]);
				//cvSetReal3D(out, x, y, i, w2[x]);
			}
		}
		for(int x=in->width-1; x>=0; x--)
		{
			elem.val[0] = w2[0*(in->width+2)+x];
			elem.val[1] = w2[1*(in->width+2)+x];
			elem.val[2] = w2[2*(in->width+2)+x];
			cvSet2D(out, x, y, elem);
		}
	}
}

void rescale_range(CvMat* img)
{
	CvScalar mean, std_dev, max_val, min_val, range;
	cvAvgSdv(img, &mean, &std_dev);

	for(int i=0; i<3; i++)
	{
		max_val.val[i] = mean.val[i] + 1.2 * std_dev.val[i];
		min_val.val[i] = mean.val[i] - 1.2 * std_dev.val[i];
		range.val[i] = max_val.val[i] - min_val.val[i];

		if(!range.val[i]) range.val[i] = 1.0;
		range.val[i] = 255./range.val[i];
	}
	
	CvMat *zeroMat = cvCreateMat(img->height, img->width, CV_64FC3);
	cvSetZero(zeroMat);
	
	cvSubS(img, min_val, img);
	cvScaleAdd(img, range, zeroMat, img);
}


int main()
{	
	char* imgName = new char[100];
	scanf("%s", imgName);
	clock_t start;
	clock_t end;
	
	IplImage* in = cvLoadImage(imgName, 1);
	//IplImage* ins[3];
	IplImage* in1 = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 3);
	cvConvertScale(in, in1, 1.);
	//ins[0] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//ins[1] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//ins[2] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	CvMat* in_mat = cvCreateMat(in->height, in->width, CV_64FC3);
	//IplImage* dst = cvCreateImage(cvSize(in->height,in->width), IPL_DEPTH_64F, 1);
	CvMat* dst = cvCreateMat(in->width,in->height, CV_64FC3);
	//IplImage* msr = cvCreateImage(cvSize(in->width,in->height),IPL_DEPTH_64F,3);
	CvMat* out = cvCreateMat(in->height,in->width, CV_64FC3);
	CvMat* msr = cvCreateMat(in->height,in->width, CV_64FC3);
	cvZero(msr);
	//IplImage* msrs[3];
	//msrs[0] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//msrs[1] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//msrs[2] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//cvZero(msrs[0]);
	//cvZero(msrs[1]);
	//cvZero(msrs[2]);

	//IplImage* outs[3];
	//outs[0] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//outs[1] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	//outs[2] = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_64F, 1);
	const int nFilter = 3;
	double sigma[nFilter];
	int default_scale = 240;
	decide_scale(sigma, nFilter, default_scale);
	//recursive gaussian
	double t=0.0;
	for(int i=0; i<nFilter; i++)
	{
		cvConvert(in1, in_mat);
		double b[5];
		//cvSplit(in1, ins[0], ins[1], ins[2], NULL);
		generate_coef3(b, sigma[i]);
		cvZero(dst);
		start = clock();
		gaussian_conv(b, in_mat, dst);
		gaussian_conv(b, dst, out);
		end = clock();
		t += double(end-start);
		printf("%lf ms elapsed\n", double(end-start));
		cvConvertScale(in_mat, in_mat, 1., 1.);
		cvConvertScale(out, out, 1., 1.);
		cvLog(in_mat, in_mat);
		cvLog(out, out);
		cvSub(in_mat, out, out);
		cvConvertScale(out, out, 1/3.);
		cvAdd(out, msr, msr);
		/*
		for(int j=0; j<3; j++) {
			cvSetZero(dst);
			start = clock();
			gaussian_conv(b, ins[j], dst);
			gaussian_conv(b, dst, outs[j]);
			end = clock();
			t += double(end-start);
			printf("%lf ms elapsed \n", double(end-start));
			cvConvertScale(ins[j],ins[j],1.,1);
			cvConvertScale(outs[j],outs[j],1.,1);
			cvLog(ins[j],ins[j]);
			cvLog(outs[j],outs[j]);
			cvSub(ins[j],outs[j],outs[j]);
			cvConvertScale(outs[j],outs[j],1/3.);
			cvAdd(outs[j],msrs[j], msrs[j]);
		}
		*/
	}
	printf("total %lf ms elapsed\n", t);
	//printf("%lf ms elapsed \n", t);

	//cvMerge(msrs[0],msrs[1],msrs[2],NULL,msr);
	rescale_range(msr);
	IplImage* msr1 = cvCreateImage(cvSize(in->width,in->height), IPL_DEPTH_8U, 3);
	cvConvertScale(msr, msr1, 1);
	cvNamedWindow("in", 1);
	cvNamedWindow("out", 1);
	cvShowImage("in", in);
	cvShowImage("out", msr1);

	
	//printf("%lf ms elapsed\n", double(end-start));
	cvWaitKey(0);
	cvDestroyAllWindows();
	

	return 0;
}